

```{r}

rm(list = ls())
library(pacman)
p_load("readr", "raster", "readxl", "dplyr", "tidyr", "rgdal", "maptools", "ggplot2", "ggthemes", "rgeos", "sp", "stringr","tidyverse", "plm", "hrbrthemes", "viridis", "forcats", "lattice", "stargazer", "dotwhisker", "lmtest", "cowplot", "ggpubr", "plyr")


Panel_long <- read_xlsx("C:/Users/tanay/Dropbox/ADS/Panel_28Jan2024.xlsx") 

```


##Fig 2
```{r}

ADS_total <-  dplyr::group_by(Panel_long, year) %>%  dplyr::summarise(ADS_total = sum(ADS2)) %>% ungroup()


extrafont::loadfonts(device="win") 

ADS_tot_plot <-  ADS_total %>% ggplot( aes(x=year, y=ADS_total)) +
    geom_line() + ylab("No. of active ADS") + xlab("Year") +
    geom_point(size = 5) + theme_clean(base_size = 16) + theme(axis.text.x=element_text(size=18), axis.text.y=element_text(size=16), axis.title.x = element_text(size=18), axis.title.y = element_text(size=16))


ggsave("C:/Users/tanay/Dropbox/ADS/Plots/Active ADS by year.png", scale = 5, ADS_tot_plot)


```


##Fig A4
```{r}
Panel_tmp <- group_by(Panel_long, year) %>% summarise(
                                    Mean_avgNL = mean(meanNL, na.rm = T), 
                                    Decile_first_avgNL = quantile(meanNL, probs = 0.1, na.rm = T),                                          Decile_last_avgNL = quantile(meanNL, probs = 0.9, na.rm = T),
                                
                                    Mean_maxNL = mean(maxNL, na.rm = T), 
                                    Decile_first_maxNL = quantile(maxNL, probs = 0.1, na.rm = T),                                           Decile_last_maxNL = quantile(maxNL, probs = 0.9, na.rm = T),
                                    
                                    
                                    Mean_dist = mean(dist, na.rm = T), 
                                    Decile_first_dist = quantile(dist, probs = 0.1, na.rm = T),                                             Decile_last_dist = quantile(dist, probs = 0.9, na.rm = T),
                                    
                                    
                                    Mean_Pop = mean(Pop, na.rm = T), 
                                    Decile_first_Pop = quantile(Pop, probs = 0.1, na.rm = T),                                               Decile_last_Pop = quantile(Pop, probs = 0.9, na.rm = T),
                                                   
                                    
                                    Mean_prop.Eleph250 = mean(prop.Eleph250, na.rm = T), 
                                    Decile_first_prop.Eleph250 = quantile(prop.Eleph250, probs = 0.1, na.rm = T),                                          
                                    Decile_last_prop.Eleph250 = quantile(prop.Eleph250, probs = 0.9, na.rm = T),
                                    
                                    
                                    Mean_prop.Agri250 = mean(prop.Agri250, na.rm = T), 
                                    Decile_first_prop.Agri250 = quantile(prop.Agri250, probs = 0.1, na.rm = T),                                          
                                    Decile_last_prop.Agri250 = quantile(prop.Agri250, probs = 0.9, na.rm = T))


NL_plot <- Panel_tmp %>% ggplot(aes(x = year)) + geom_line(aes(y = Mean_avgNL), show.legend = T) + geom_point(aes(y = Mean_avgNL, shape = "Mean_avgNL"), size = 2) + geom_line(aes(y = Decile_first_avgNL), show.legend = T) + geom_point(aes(y = Decile_first_avgNL, shape = "Decile_first_avgNL"), size = 2) + geom_line(aes(y = Decile_last_avgNL), show.legend = T) + geom_point(aes(y = Decile_last_avgNL, shape = "Decile_last_avgNL"), size = 2) +  scale_shape_manual("", labels = c("Mean", "10th %tile", "90th %tile"), breaks = c("Mean_avgNL", "Decile_first_avgNL", "Decile_last_avgNL"), values = c(16, 17, 18)) + labs(y = "Mean_NL (cd)", x = "Year") + theme_clean(base_size = 20)

Dist_plot <- Panel_tmp %>% ggplot(aes(x = year)) + geom_line(aes(y = Mean_dist), show.legend = T) + geom_point(aes(y = Mean_dist, shape = "Mean_dist"), size = 2) + geom_line(aes(y = Decile_first_dist), show.legend = T) + geom_point(aes(y = Decile_first_dist, shape = "Decile_first_dist"), size = 2) + geom_line(aes(y = Decile_last_dist), show.legend = T) + geom_point(aes(y = Decile_last_dist, shape ="Decile_last_dist"), size = 2) + scale_shape_manual("", labels = c("Mean", "10th %tile", "90th %tile"), breaks = c("Mean_dist", "Decile_first_dist", "Decile_last_dist"), values = c(16, 17, 18)) +  labs(y = "Dist_eleph_habitat (m)", x = "Year") + theme_clean(base_size = 20)

Pop_plot <- Panel_tmp %>% ggplot(aes(x = year)) + geom_line(aes(y = Mean_Pop), show.legend = T) + geom_point(aes(y = Mean_Pop, shape = "Mean_Pop"), size = 2) + geom_line(aes(y = Decile_first_Pop), show.legend = T) + geom_point(aes(y = Decile_first_Pop, shape = "Decile_first_Pop"), size = 2) + geom_line(aes(y = Decile_last_Pop), show.legend = T) + geom_point(aes(y = Decile_last_Pop, shape = "Decile_last_Pop"), size = 2) + scale_shape_manual("", labels = c("Mean", "10th %tile", "90th %tile"), breaks = c("Mean_Pop", "Decile_first_Pop", "Decile_last_Pop"), values = c(16, 17, 18)) +  labs(y = "Population (count)", x = "Year") + theme_clean(base_size = 20)

prop.Eleph250_plot <- Panel_tmp %>% ggplot(aes(x = year)) + geom_line(aes(y = Mean_prop.Eleph250), show.legend = T) + geom_point(aes(y = Mean_prop.Eleph250, shape = "Mean_prop.Eleph250"), size = 2) + geom_line(aes(y = Decile_first_prop.Eleph250), show.legend = T) + geom_point(aes(y = Decile_first_prop.Eleph250, shape = "Decile_first_prop.Eleph250"), size = 2) + geom_line(aes(y = Decile_last_prop.Eleph250), show.legend = T) + geom_point(aes(y = Decile_last_prop.Eleph250, shape ="Decile_last_prop.Eleph250"), size = 2) + scale_shape_manual("", labels = c("Mean", "10th %tile", "90th %tile"), breaks = c("Mean_prop.Eleph250", "Decile_first_prop.Eleph250", "Decile_last_prop.Eleph250"), values = c(16, 17, 18)) +  labs(y = "Prop_eleph_250", x = "Year") + theme_clean(base_size = 20)


Stats <- ggarrange(NL_plot, Dist_plot, Pop_plot, prop.Eleph250_plot, nrow = 2, ncol = 2, common.legend = T)


```



##Fig A5
```{r}

Panel_tmp <-  Panel_long %>%
  group_by(year) %>%
  mutate(HDI_total = sum(HDI)) %>% ungroup()

HDI_total <- dplyr::select(Panel_tmp, year, HDI_total)
HDI_total <- HDI_total[!duplicated(HDI_total), ]

HDI_tot_plot <- HDI_total %>%
  ggplot( aes(x=year, y=HDI_total)) +
    geom_line() +
    geom_point(size=3) + ylab("No. of human deaths") + xlab("Year") +
    theme_clean(base_size = 18) + theme(axis.text.x=element_text(size=16), axis.text.y=element_text(size=16), legend.text = element_text(size = 14))

Panel_tmp <-  Panel_long %>%  mutate(HDI_RK = ifelse(!is.na(HDI_0_1), ifelse(EDI> 0 & HDI_0_1 > 0, HDI_0_1, 0), NA)) %>% ungroup() %>% group_by(year) %>%
  mutate(EDI_total = sum(EDI), HD_0_1_sum = sum(HDI_RK), EDI_accident = sum(EDI_accident), EDI_other = sum(EDI_other)) %>% ungroup()

EDI_total <- dplyr::select(Panel_tmp, year, EDI_total, HD_0_1_sum, EDI_accident, EDI_other)

EDI_total <- EDI_total[!duplicated(EDI_total), ]

EDI_tot = pivot_longer(EDI_total, cols = c("EDI_total", "EDI_accident", "EDI_other", "HD_0_1_sum"), names_to = "type", values_to = "EDI")

EDI_tot$type <- gsub(c("EDI_accident", "EDI_other"), c("EDI_PA", "EDI_O"), EDI_tot$type)

EDI_decomp_plot <- ggplot(data = EDI_tot, aes(x = year, y = EDI, group = type)) + ylab("No. of elephant deaths") + xlab("Year") + geom_line() + geom_point(aes(shape = type), size = 4) + scale_shape_manual(name = "", values = c(15, 19, 17, 18)) + theme_clean(base_size = 18) + theme(legend.key.size = unit(0.4, "mm"), legend.position = c(0.65, 0.75), axis.text.x=element_text(size=16), axis.text.y=element_text(size=16), legend.text = element_text(size = 14), legend.background = element_rect(color = NA))


HD_ED <- ggarrange(HDI_tot_plot, EDI_decomp_plot, nrow = 2, ncol = 1)


```

